{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Temperature-Dependent Power Flow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement Temperature-Dependent Power Flow (TDPF) in pandapower based on the following publications:\n",
    "\n",
    "S. Frank, J. Sexauer and S. Mohagheghi, \"Temperature-dependent power flow\", IEEE Transactions on Power Systems, vol. 28, no. 4, pp. 4007-4018, Nov 2013.\n",
    "\n",
    "B. Ngoko, H. Sugihara and T. Funaki, \"A Temperature Dependent Power Flow Model Considering Overhead Transmission Line Conductor Thermal Inertia Characteristics,\" 2019 IEEE International Conference on Environment and Electrical Engineering and 2019 IEEE Industrial and Commercial Power Systems Europe (EEEIC / I&CPS Europe), 2019, pp. 1-6, doi: 10.1109/EEEIC.2019.8783234.\n",
    "\n",
    "This tutorial demonstrates how to use this feature in pandapower in order to calculate temperature of overhead lines and also obtain the results of the TDPF with the adjusted line parameters. Furthermore, thermal inertia of overhead lines can be considered by specifying the time delay after which the temperature should be calculated. This implementation also supports the distributed slack power flow calculation.\n",
    "\n",
    "We use the approach from Frank et al. to calculate the thermal resistance of the lines and to define the Jacobian matrix and the mismatch eqiation. The approach of Ngoko et al. is used to update the thermal resistance based on the line current and the weather parameters, and to introduce the thermal delay term in the Jacobian matrix. This allows using both methods. Furthermore, the approach to use thermal resistance is more versatile because it allows expanding the implementation for underground cables and transformers in the future.\n",
    "\n",
    "This implementation will be of benefit for grid studies of Dynamic Line Rating (DLR). In addition, a calculation of a grid state at a specified time after an N-1 event provides additional capabilities for contingency studies and for obtaining emergency line ratings for a given grid. \n",
    "\n",
    "In this tutorial, we demonstrate the use of this feature with a small example.\n",
    "\n",
    "After importing the necessary libraries, we create a 5-bus power system (inspired by Ngoko et al.) that we use to demonstrate the calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:26.149166Z",
     "start_time": "2025-10-20T13:00:07.256475Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pandapower.test.loadflow.test_tdpf import simple_test_grid\n",
    "from pandapower.pf.create_jacobian_tdpf import calc_T_ngoko, calc_a0_a1_a2_tau, calc_T_frank, calc_g_b, calc_i_square_p_loss, calc_r_theta, calc_r_theta_from_t_rise\n",
    "from pandapower.pypower.idx_brch import BR_R, BR_X\n",
    "from pandapower.plotting.simple_plot import simple_plot\n",
    "from pandapower.run import runpp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:28.130896Z",
     "start_time": "2025-10-20T13:00:26.183858Z"
    }
   },
   "outputs": [],
   "source": [
    "net = simple_test_grid(load_scaling=0.25, sgen_scaling=0.5, with_gen=True, distributed_slack=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:28.797855Z",
     "start_time": "2025-10-20T13:00:28.163710Z"
    }
   },
   "outputs": [],
   "source": [
    "simple_plot(net, plot_loads=True, plot_sgens=True, plot_gens=True, load_size=3, sgen_size=3, gen_size=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:38.706757Z",
     "start_time": "2025-10-20T13:00:28.851507Z"
    }
   },
   "outputs": [],
   "source": [
    "runpp(net)\n",
    "net.res_line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:39.642550Z",
     "start_time": "2025-10-20T13:00:39.600194Z"
    }
   },
   "outputs": [],
   "source": [
    "net.res_bus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculation of the overhead line temperature\n",
    "\n",
    "It is possible to use the functions from the implementation separately, in order to calculate the overhead line temperature.\n",
    "\n",
    "The implementation in Ngoko et al. defines three terms to represent a simplified thermal model of an overhead line. The parameters $a_0$, $a_1$, and $a_2$ specify a constant, linear, and quadratic terms of the dependence of the line temperature from the square of the current $I$ that is flowing through the overhead line:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "T_{ss} = a_0 + a_1 \\cdot I^2 + a_2 \\cdot I^4\n",
    "$$\n",
    "\n",
    "$$\n",
    "T = T_{ss} - (T_{ss} - T_0) \\cdot \\exp(-t/\\tau)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the equations above, $T_{ss}$ stands for the steady-state temperature, $T_0$ is the initial temperature, and $T$ stands for the temperature at the time $t$ after a step change of current. The parameter $\\tau$ is a time constant, used in the caluclation of line temperature with thermal inertia. After $\\tau$ seconds, the temperature reaches approximately 63.2 % of the steady-state value.\n",
    "\n",
    "This approach requires the following inputs (default values are in parentheses):\n",
    "\n",
    "* air temperature (20 °C)\n",
    "* initial line temperature (20 °C)\n",
    "* reference temperature (20 °C)\n",
    "* resistance at reference temperature (r_ohm_per_km)\n",
    "* conductor outer diameter\n",
    "* specific mass of the conductor multiplied by the specific thermal capacity of the material\n",
    "* wind speed (0.6 m/s)\n",
    "* wind angle of attack (45 °)\n",
    "* solar radiation (900 W/m²)\n",
    "* solar emissivity and solar absorptivity (0.5)\n",
    "* thermal coefficient of resistivity (4.03e-3)\n",
    "\n",
    "The parameters are specified in the net.line element table. If any parameters are missing, they will be substitutet by default assumptions. The parameter \"tdpf\" is always required. If the option to not update $R_\\Theta$ is used, then \"r_theta_kelvin_per_mw\" is required. Otherwise, \"conductor_outer_diameter_m\" is required. Furthermore, if thermal inertia is considered by setting the option \"tdpf_delay_s\", the parameter \"mc_joule_per_m_k\" is required."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:40.814662Z",
     "start_time": "2025-10-20T13:00:40.780487Z"
    }
   },
   "outputs": [],
   "source": [
    "t_air_pu = 35\n",
    "alpha_pu = 4.03e-3\n",
    "r_ref = net.line.r_ohm_per_km.values / 1e3\n",
    "a0, a1, a2, tau = calc_a0_a1_a2_tau(t_air_pu, 80, 20, r_ref, 30.6e-3,\n",
    "                                    1490, 0.6, 45, 900, alpha_pu, 0.5, 0.5)\n",
    "calc_T_ngoko(np.square(net.res_line.i_ka.values * 1e3), a0, a1, a2, None, None, None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the approach of Frank et al., the parameters specifying the thermal model according to Ngoko et al. are used to calculate the thermal resistamce $R_{\\Theta}$: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "R_{\\Theta} = \\frac{T_{Rise}}{P_{Loss}} = \\frac{a_0 - T_{air} + a_1 \\cdot I^2 + a_2 \\cdot I^4}{P_{Loss}}\n",
    "$$\n",
    "\n",
    "$$\n",
    "T = T_{air} + R_\\Theta \\cdot P_{Loss}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, the thermal resistance can be approximated using an assumption for the rated temperature rise, as described in Frank et al. To this end, the air temperature $T_{air}$ and the temperature rise are added to calculate the line temperature. The temperature rise is calculated by multiplying the thermal resistance $R_\\Theta$ and the active power losses $P_{Loss}$. The thermal resistance is included in the line results table."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:41.096161Z",
     "start_time": "2025-10-20T13:00:41.065315Z"
    }
   },
   "outputs": [],
   "source": [
    "branch = net._ppc[\"branch\"]\n",
    "tdpf_lines = np.ones(len(branch)).astype(bool)\n",
    "r = branch[tdpf_lines, BR_R].real\n",
    "x = branch[tdpf_lines, BR_X].real\n",
    "g, b = calc_g_b(r, x)\n",
    "Vm = abs(net._ppc[\"internal\"][\"V\"])\n",
    "Va = np.angle(net._ppc[\"internal\"][\"V\"])\n",
    "i_square_pu, p_loss_pu = calc_i_square_p_loss(branch, tdpf_lines, g, b, Vm, Va)\n",
    "r_theta_pu = calc_r_theta(t_air_pu, a0, a1, a2, np.square(net.res_line.i_ka.values*1e3), p_loss_pu)\n",
    "calc_T_frank(p_loss_pu, t_air_pu, r_theta_pu, None, None, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:41.206728Z",
     "start_time": "2025-10-20T13:00:41.171501Z"
    }
   },
   "outputs": [],
   "source": [
    "# using an approximation for the rated temperature rise:\n",
    "r_theta_kelvin_per_mw = calc_r_theta_from_t_rise(net, 25)\n",
    "calc_T_frank(p_loss_pu, t_air_pu, r_theta_kelvin_per_mw.values * net.sn_mva / 1, None, None, None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By providing the option tdpf=True, the overhead line temperature can be obtained from the power flow caluclation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:48.614953Z",
     "start_time": "2025-10-20T13:00:41.260766Z"
    }
   },
   "outputs": [],
   "source": [
    "runpp(net, tdpf=True, max_iteration=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taking the results of the power flow calculation for the line current and using them as the input for the introduced functions, we can make sure that the temperature calculation delivers the same results. However, if we repeat the power flow calculation without the TDPF, the results will be different because the effect of the temperature on the resitance and the current is ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:48.721158Z",
     "start_time": "2025-10-20T13:00:48.648916Z"
    }
   },
   "outputs": [],
   "source": [
    "net.res_line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:49.060721Z",
     "start_time": "2025-10-20T13:00:49.037031Z"
    }
   },
   "outputs": [],
   "source": [
    "calc_T_ngoko(np.square(net.res_line.i_ka.values * 1e3), a0, a1, a2, None, None, None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:49.727720Z",
     "start_time": "2025-10-20T13:00:49.524843Z"
    }
   },
   "outputs": [],
   "source": [
    "runpp(net)\n",
    "calc_T_ngoko(np.square(net.res_line.i_ka.values * 1e3), a0, a1, a2, None, None, None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculation works for distributed slack power flow, too. However, the parameter net.sn_mva should be increased in order to avoid numerical instability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:54.650295Z",
     "start_time": "2025-10-20T13:00:50.118121Z"
    }
   },
   "outputs": [],
   "source": [
    "net.sn_mva = 1000\n",
    "runpp(net, tdpf=True, tdpf_delay_s=5 * 60, distributed_slack=True, max_iteration=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we demonstrate the calculation with a time delay of 5 minutes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:55.150475Z",
     "start_time": "2025-10-20T13:00:54.674366Z"
    }
   },
   "outputs": [],
   "source": [
    "runpp(net, tdpf=True, tdpf_delay_s=5 * 60, max_iteration=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:00:55.270468Z",
     "start_time": "2025-10-20T13:00:55.175412Z"
    }
   },
   "outputs": [],
   "source": [
    "net.res_line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the temperature rise is substantially lower due to thermal inertia. The thermal inertia effect can be illustrated in the following figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-10-20T13:01:09.089728Z",
     "start_time": "2025-10-20T13:00:55.614187Z"
    }
   },
   "outputs": [],
   "source": [
    "delays = np.arange(0, 65, 5)\n",
    "delay_tab = pd.DataFrame(index=delays, columns=net.line.index.values)\n",
    "\n",
    "for d in delays:\n",
    "    runpp(net, tdpf=True, tdpf_delay_s=d * 60, max_iteration=30)\n",
    "    delay_tab.loc[d, :] = net.res_line.temperature_degree_celsius.values\n",
    "    \n",
    "delay_tab.plot(ylabel=\"Temperature (°C)\", xlabel=\"Time delay (min)\", \n",
    "               title=\"Time delay and overhead line temperature\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
